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Summary 


The  goal  of  this  project  was  to  develop  new  mathematical  and  computational  techniques  for 
quantifying  the  errors  in  seismic  event  locations,  with  the  goal  of  improving  on  the  epicentral 
confidence  ellipses  and  focal  depth  confidence  intervals  in  conventional  use.  The  primary  focus 
of  the  project  was  to  develop  a  rigorous  approach  to  the  effects  of  model  errors,  or  errors  in 
the  travel-time  predictions  calculated  from  an  assumed  velocity  model  for  the  Earth.  The 
approach  that  was  investigated  associates  model  errors  with  the  uncertainty  in  path  travel¬ 
time  corrections  that  are  inferred  from  a  calibration  analysis.  This  naturally  leads  to  the  idea 
of  analyzing  location  uncertainty  in  the  context  of  a  joint  inverse  problem  that  treats  event 
location  and  travel-time  calibration  as  coupled  problems. 

A  primary  accomplishment  of  the  project  was  a  general  theoretical  formulation  of  the  joint 
location/calibration  inverse  problem  that  applies  to  a  wide  range  of  situations,  including  the  use 
of  non-Gaussian  distributions  for  pick  errors,  the  accommodation  of  non-GTO  calibration  events, 
and  the  use  of  a  general  parameterization  of  travel-time  corrections  that  subsumes  a  variety  of 
calibration  techniques  previous  researchers  have  considered.  The  latter  includes  simple  station- 
specific  time  terms,  station-specific  travel-time  correction  surfaces,  and  3D  velocity  models. 
The  formulation  leads  to  a  general  expression  for  event  location  confidence  regions  based  on 
the  statistical  framework  of  hypothesis  testing  with  likelihood-ratio  statistics. 

A  second  accomplishment  of  the  project  was  to  implement  the  theoretical  formulation  for 
the  case  of  time-term  corrections  as  a  proof  of  concept  of  the  joint  inversion  uncertainty  ap¬ 
proach.  Applications  performed  with  data  from  Nevada  Test  Site  explosions  illustrate  the  effects 
of  model  errors  on  event  location  confidence  regions  and  their  dependence  on  the  location  ac¬ 
curacy  of  calibration  events.  Despite  attempts  to  increase  the  computational  efficiency  of  the 
joint  inversion  approach,  these  applications  also  indicated  that  the  approach  is  not  feasible  for 
complex  methods  of  travel-time  calibration  like  3-D  velocity  tomography. 

The  third  accomplishment  of  the  project  was  the  development  of  an  approximate,  but  much 
more  efficient,  version  of  the  new  uncertainty  approach  that  follows  the  traditional  factoriza¬ 
tion  of  travel-time  calibration  and  event  location  into  a  two-stage  process.  The  approximation 
assumes  that  the  uncertainty  in  travel-time  corrections  estimated  in  a  calibration  analysis  can 
be  represented  by  a  Gaussian  probability  distribution.  The  generality  of  the  theoretical  formu¬ 
lation  is  otherwise  preserved.  To  test  the  two-stage  approach,  new  algorithms  were  developed 
for  calculating  calibration  uncertainty  (stage  one)  and  employing  it  in  the  calculation  of  loca¬ 
tion  confidence  regions  (stage  two).  Application  of  these  algorithms  to  the  Nevada  Test  Site 
data,  assuming  time-term  corrections,  demonstrate  the  much  greater  efficiency  of  the  two-stage 
approach  and  suggest  that  the  approach  is  feasible  for  the  case  of  tomographic  calibration. 
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Chapter  1 

Introduction 


The  assessment  of  errors  in  seismic  event  locations  remains  a  challenging  problem  for  nuclear 
monitoring  research.  Theoretical  treatments  of  location  uncertainly  have  generally  focused 
on  the  effects  of  observational,  or  pick,  errors  in  the  seismic  arrival  times  (or  other  observed 
attributes)  used  in  locating  an  event.  However,  it  has  been  recognized  for  a  long  time  that 
an  equal  or  larger  source  of  location  uncertainty,  especially  when  arrivals  at  regional  distances 
are  involved,  is  the  error  in  the  travel-time  predictions  that  are  calculated  with  models  of  the 
Earth’s  velocity  structure.  Formulations  of  location  uncertainty  have  generally  lacked  a  rigorous 
treatment  of  these  model  errors,  using  instead  the  ad  hoc  assumption  that  they  simply  augment 
the  pick  errors. 

The  project  reported  here  has  attempted  to  develop  a  rigorous,  and  general,  treatment  of 
model  errors  by  linking  them  to  the  errors  incurred  in  the  calibration  of  seismic  travel-times. 
This  association  leads  naturally  to  consideration  of  a  joint  inverse  problem  in  which  arrival 
time  data  from  multiple  events  are  used  to  infer  the  event  locations  together  with  calibration 
parameters  that  generate  corrections  to  the  model-based  travel-time  predictions.  Taking  one 
of  the  multiple  events  to  be  the  target  of  investigation,  and  the  others  as  historical  calibration 
events,  an  error  analysis  performed  on  the  joint  location/calibration  inverse  problem  implicitly 
accounts  for  the  effects  of  model  errors  on  the  target  event  location.  The  joint  inversion  approach 
also  provides  a  means  to  address  the  important  issue  of  how  errors  in  the  locations  of  calibration 
events  influence  the  location  uncertainty  for  the  target  event. 

The  next  chapter  of  this  report  formulates  the  problem  of  event  location  uncertainty  and 
reviews  the  classical  theory  of  elliptical  confidence  regions  on  locations,  which  follows  from  the 
linearization  of  the  travel-time  forward  problem  and  the  assumption  that  pick  errors  obey  a 
Gaussian  probability  distribution.  The  linear /Gaussian  theory  is  then  generalized  based  on  the 
statistical  framework  of  maximum-likelihood  estimation  and  hypothesis  testing,  and  numerical 
techniques  for  the  evaluation  of  non-elliptical  confidence  regions  are  described. 

Chapter  3  of  the  report  extends  the  theory  and  algorithms  from  single-event  location  to 
multiple-event  location,  i.e.  the  joint  location/calibration  problem.  The  resulting  approach  is 
illustrated  with  examples  that  use  data  from  well-located  explosions  at  the  Nevada  Test  Site. 

Chapter  4  presents  a  two-stage  method  that,  following  conventional  practice,  separates  the 
calibration  and  location  aspects  of  the  joint  inverse  problem.  The  separation  is  developed 
under  a  specific  approximation  to  the  full  theory;  namely,  that  the  uncertainty  in  calibration 
parameters  can  be  represented  by  a  multivariate  Gaussian  distribution.  Appropriate  algorithms 
to  implement  the  two-stage  approach  are  presented  and  illustrated  with  the  Nevada  Test  Site 
data.  The  report  ends  with  some  conclusions  and  suggestions  for  future  work. 
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Chapter  2 


Uncertainty  Analysis  in  Single-Event 
Location 

2.1  The  Single-Event  Location  Problem  and  Approaches  to  Un¬ 
certainty 

Most  event  location  algorithms  infer  the  hypocenter  and  origin  time  of  an  event  from  a  set 
of  arrival  times  of  various  seismic  phases  observed  at  various  stations  of  a  network.  We  will 
denote  the  arrival-time  data  d%  for  1  =  1.  . . . ,  n,  where  n  is  the  total  number  of  station/phase 
combinations  that  have  been  observed.  Letting  x  denote  the  event  hypocenter  and  t  its  origin 
time,  we  can  express  the  single-event  location  problem  as 

di  =  Tj(x)  +t  +  ei:  i  =  1, . . .  ,n,  (2.1) 

where  e*  is  the  observational  error  in  di,  and  Tr  is  a  function  that  predicts  the  travel-time  of 
the  ith  station/phase  pair  as  a  function  of  the  event  hypocenter.  The  travel-time  prediction 
functions  are  usually  based  on  a  velocity  model  for  the  Earth.  In  matrix/vector  notation  we 
can  write  equation  (2.1)  as 

d  =  T(x)  +  It  +  e,  (2.2) 

where  each  element  of  the  column  vector  1  is  unity. 

A  solution  of  the  single-event  location  problem  is  commonly  taken  to  comprise  point  esti¬ 
mates  of  the  event  location  parameters,  which  we  will  denote  (x,  t),  together  with  a  probabilistic 
description  of  the  uncertainty  in  these  estimates.  The  uncertainty  description  follows  from  an 
assumed  probabilistic  model  for  the  observational  errors,  e*.  There  are  a  number  of  accepted 
ways  to  describe  uncertainty.  For  example,  if  x  is  found  to  be  a  Gaussian  random  variable,  its 
uncertainty  is  well-characterized  by  a  variance-covariance  matrix,  V,  such  that 

E[(x  -  x)(x  —  x)T]  =  V,  (2.3) 

where  E[  ]  denotes  the  mathematical  expectation  of  a  random  variable.  The  presumption  here 
is  that  x  is  an  unbiased  estimate  of  x,  or  E[x]  =  x.  There  is  no  guarantee,  however,  that  an 
unbiased  Gaussian  estimate  of  x  exists  or,  if  it  does,  that  its  variance  is  independent  of  x. 

A  more  general  description  of  uncertainty  is  accomplished  with  Neyman-Pearson  confidence 
regions  on  the  unknown  parameters  or  subsets  thereof.  A  confidence  region  on  x,  for  example, 
is  a  region  of  x-space,  calculated  from  the  observed  data,  to  which  we  can  assign  a  probability  of 
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including  the  true  hypocenter  x.  Letting  (3  be  that  probability,  or  confidence  level  (e.g.  f3  =  0.9 
for  a  90%  confidence  region),  and  letting  Xp(d)  be  its  associated  confidence  region,  we  have 

Prob  [x  e  Xp{d)]  =  p.  (2.4) 

The  probability  operator  used  here  is  that  which  is  induced  by  the  assumed  error  model  for  e 
and,  hence,  d. 

Bayesian  inference  provides  another  approach  to  uncertainty.  In  this  approach,  the  location 
parameters  are  considered  to  be  random  variables,  and  a  complete  solution  of  the  inverse 
problem  is  embodied  in  their  joint  posterior  probability  density  function  (p.d.f.),  denoted  /[x,  t  \ 
d]  (where  “|d”  indicates  conditioning  on  the  data  vector).  From  the  posterior  distribution,  a 
point  estimate  of  x  is  available  as  the  MAP  (maximum  a  posteriori)  estimate: 

/[x  |  d]  =  max/[x  |  d],  (2.5) 

X 

where  /[x  |  d]  is  the  marginal  posterior  p.d.f.  of  x,  defined  as 

/[x  |  d ]  =  J  dt  /[x,t  |  d],  (2.6) 

A  Bayesian  confidence  region  can  be  defined  from  the  posterior  p.d.f.  with 

[  dx  /[x  |  d]  =  P  (2.7) 

J  xeA’g(d) 

but  such  confidence  regions  do  not  necessarily  obey  the  confidence  property  stated  in  equation 
(2-4)- 

This  project  has  focused  on  Neyman-Pearson  confidence  regions  as  a  framework  for  event 
location  uncertainty  analysis.  Before  describing  this  approach  in  a  general  context,  we  describe 
it  for  the  special  case  most  often  assumed  in  seismic  event  location  studies. 


2.2  Gaussian/Linear  Confidence  Regions 

The  traditional  solution  to  the  single-event  location  problem  consists  of  least-squares  estimates 
of  the  event  parameters.  These  minimize  a  data-misfit  function  given  by 

n 

T(x,f;d)  =  ^Wi  (di  -  T*(x)  -  f)2  ,  (2.8) 

2—1 

where  the  Wi  are  prescribed  weights.  We  can  also  write  the  least-squares  misfit  function  as 

T(x,  t;  d)  =  (d  -  T(x)  -  It)  TW(d  -  T(x)  -  It)  (2.9) 

with  the  diagonal  matrix  W  containing  the  weights: 


The  least-squares  parameter  estimates  satisfy 

T(x,t;d)  =  min  'F(x,  f;  d). 
x.t 


(2.10) 


(2.11) 
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Numerically,  the  minimization  of  \k  is  typically  performed  with  a  Gauss-Newton  algorithm,  as 
proposed  originally  by  Geiger  (1912)  and  re-formulated  by  Flinn  (1965). 

Flinn  (1965)  also  developed  the  theory  for  location  confidence  regions  in  the  least-squares 
setting,  and  we  redevelop  this  result  here  as  a  prelude  to  a  more  general  theory.  Flinn  derives 
Neyman-Pearson  confidence  regions  in  a  hypothesis  testing  framework.  The  gist  of  this  approach 
is  as  follows.  Given  an  arbitrary  hypocenter  x,  we  take  as  a  null  hypothesis  that  the  true  event 
hypocenter  is  x.  If,  based  on  a  test  statistic  r(x,  d),  we  cannot  reject  the  null  hypothesis  at 
some  level  of  confidence,  we  take  x  to  be  a  point  inside  the  confidence  region  at  that  level. 

Flinn  applied  this  approach  with  the  underlying  assumption  that  the  observational  errors  are 
zero-mean,  Gaussian  random  variables  whose  standard  deviations  are  known  within  a  constant 
factor,  a.  The  probability  density  function  (p.d.f.)  of  e*  in  this  case  can  be  written  as 

(2-l2) 

where  i ^  is  a  nominal  standard  deviation  assigned  to  el.  With  this  assumption,  and  assum¬ 
ing  that  the  errors  are  statistically  independent  (uncorrelated)  the  weights  for  the  data-misfit 
function  are  naturally  set  as 

Wi  ~  v~2.  (2.13) 


Location  confidence  regions  can  be  defined  on  all  the  event  location  parameters  (x  and  t)  in 
four-dimensional  space,  or  on  any  subset  of  the  parameters.  Here  we  will  consider  confidence 
regions  on  the  hypocenter  x,  thus  treating  t  as  a  “nuisance”  parameter.  For  these,  Flinn  (1965) 
applied  hypothesis  testing  with  the  test  statistic  given  by 


r(x,d) 


—  min  T(x,  t;  d)  —  \H(x,  t;  d)  , 


(2.14) 


where  (x,  t)  is  the  least-squares  location  solution  and  a  is  an  estimate  of  a  which  Flinn  took  to 
be 


a2  = 


1 


n  —  4 


T(£,?;d). 


(2.15) 


We  see  that  r  compares  the  data  misfit  achieved  by  a  given  hypocenter  x  (allowing  t  to  adjust) 
to  the  minimum  misfit  achieved  by  any  x.  A  confidence  region  on  x  at  the  /3  confidence  level, 
or  Ag(d),  comprises  those  points  x  yielding  values  of  the  test  statistic  below  some  cutoff.  That 
is, 


Xp{d)  =  {x:  t(x,  d)  <  rg}  , 


(2.16) 


where  Tg  is  a  critical  value  of  the  test  statistic  satisfying 


Prob 


r(x,  d)  <  Tg 


(3- 


(2.17) 


Equation  (2.16)  can  be  interpreted  to  mean  that  x  can  be  rejected  as  the  true  hypocenter  if 
its  relative  data  misfit,  as  measured  by  r(x,  d),  is  too  large.  Equation  (2.17)  states  that  the 
cutoff  value  for  this  rejection,  Tg,  marks  the  /3  point  on  the  cumulative  probability  distribution 
of  r(x,  d),  as  determined  under  the  null  hypothesis.  This  choice  for  the  cutoff  ensures  the 
confidence  property  of  equation  (2.4).  That  is,  the  true  hypocenter  will  be  rejected  with  prob¬ 
ability  1  —  (3  or,  stated  conversely,  the  confidence  region  will  include  the  true  hypocenter  with 
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probability  (3.  This  interpretation  presumes,  however,  that  the  distribution  of  r(x,  d)  does  not 
depend  on  the  parameters  t  and  a,  which  we  will  see  is  indeed  the  case  in  the  Gaussian/linear 
analysis. 

To  obtain  the  confidence  region  on  x  in  a  convenient  form  we  need  a  numerical  value  for 
T/3  and  a  way  to  parameterize  the  points  satisfying  (2.16).  These  requirements  are  easily  met 
under  the  Gaussian  error  assumption  if  we  make  the  additional  assumption  that  the  travel-time 
functions  can  be  well-approximated  by  a  linear  expansion  about  the  least-squares  estimate: 

T(x)  =  T(x)  +  A(x  —  x).  (2.18) 


Here,  A  is  a  Jacobian  matrix  containing  derivatives  of  the  travel-time  functions,  evaluated  at 
x  =  x: 


Aij  — 


dTi  (x) 


d'ju 


(2.19) 


Given  (2.18),  the  location  problem  becomes  a  linear  inverse  problem  (the  dependence  of  t  is 
exactly  linear),  and  standard  methods  of  Gaussian/linear  analysis  can  be  applied  to  determine 
t  and  its  probability  distribution  (e.g.,  Anderson,  1965).  We  report  the  results. 

Under  the  linearization  of  T,  the  test  statistic  can  be  expressed  as 


1 


r(x,  d)  =  —  (x  —  x) 1  A 1  Q  WQA(x-x), 


(j 


where  Q  is  the  projection  matrix  given  by 
1 


Q  =  I- 


uwi 


11TW. 


(2.20) 


(2.21) 


Given  that  e  is  Gaussian,  r  is  the  (scaled)  ratio  of  two  independent  chi-squared  random  variables 
having  3  and  n  —  4  degrees  of  freedom,  respectively,  r  itself  (divided  by  3)  is  thus  F  distributed 
with  these  degrees  of  freedom,  and  we  find  that 


r/3  =  3  Fg(3,  ra  —  4).  (2.22) 

Ff3(k,m)  denotes  the  /3-point  of  the  F  distribution  for  k  and  m  degrees  of  freedom.  The 
confidence  region  on  x  at  confidence  level  /3  is  thus  given  by 

(x  —  x)TATQrWQA(x  —  x)  <  3  a2  Fp(3,n -4).  (2.23) 

We  see  from  equation  (2.22)  that  does  not  depend  on  any  of  the  unknown  parameters 
(x,  t,  a).  Equation  (2.23)  therefore  describes  the  interior  of  an  ellipsoid  in  hypocenter  space, 
centered  at  the  least-squares  estimate  x.  The  parameters  of  the  ellipsoid  can  be  determined 
from  the  eigenvectors  and  eigenvalues  of  the  matrix  A1  Q1  WQA.  The  inverse  square-root  of 
the  eigenvalues,  multiplied  by  the  square-root  of  the  right-hand  side  of  (2.23),  are  the  semi-axis 
lengths  of  the  ellipsoid. 

An  alternative  to  Flinn’s  (1965)  results  was  proposed  by  Evernden  (1969)  for  the  case  in 
which  a  is  taken  to  be  a  known  quantity.  In  this  case,  t(x,  d)  in  equation  (2.20),  with  a  set 
to  the  known  a,  becomes  chi-squared  distributed  with  3  degrees  of  freedom  and  the  confidence 
region  becomes 

(x  —  x)TArQI'WQA(x  —  x)  <  ct2x|(3).  (2.24) 

The  confidence  region  scaling  is  smaller  when  the  chi-squared  distribution  is  used  for  r,  sig¬ 
nifying  that  a  loss  of  information  occurs  when  a  is  not  known  and  has  to  be  estimated  from 
the  data.  Table  2.1  shows  the  ratio  of  confidence  region  scaling  factors  implied  by  Flinn’s  and 
Evernden’s  methods. 
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Table  2.1:  Ratio  of  F  to  x 2  scaling  of  hypocenter  confidence  regions 


Ratio  = 


r  i  i/2 

3Fp(3,  n  —  4) 

Xj|(3) 


n 

(3  =  90% 

/?  =  95% 

5 

5.071 

9.100 

6 

2.097 

2.712 

8 

1.418 

1.591 

10 

1.256 

1.351 

12 

1.185 

1.249 

15 

1.130 

1.174 

20 

1.087 

1.115 

30 

1.052 

1.069 

50 

1.029 

1.038 

2.3  Non-Elliptical  Confidence  Regions 

In  previous  work  (e.g.  Rodi  and  Toksoz,  2000),  we  generalized  the  linear /Gaussian  theory 
for  location  confidence  regions  to  allow  non-Gaussian  probability  distributions  for  pick  errors, 
to  accommodate  nonlinear  constraints  on  location  parameters,  and  to  take  account  of  the 
nonlinearity  of  the  travel-time  functions.  We  summarize  this  generalization  here. 

2.3.1  Likelihood  function  for  arrival  data 

We  still  assume  the  pick  errors  are  statistically  independent  but,  following  Billings  et  al.  (1994), 
allow  each  to  be  distributed  with  a  generalized  Gaussian  distribution,  whose  probability  density 
function  (p.d.f.)  is  given  by 

/fo]  =  JFT\  exP  i~-  ~  ),  P>  L  (2-25) 

K(p)Ui  {  P  CTi  ) 

where  at  is  a  standard  error  and 

K(p)  =  2p1/pT(l  +  1  Ip).  (2.26) 

(T  is  the  gamma  function.)  When  p  =  2,  the  p.d.f.  is  Gaussian  (normal)  with  a  mean  of  zero 
and  variance  of  of.  When  p  =  1,  /[  ]  is  a  Laplace  distribution  (two-sided  exponential).  Like 
we  did  earlier,  we  assume  that  the  standard  errors  are  known  within  a  scale  parameter  a: 

Vi  =  avi  (2.27) 

with  the  being  known,  nominal  standard  errors. 

Given  that  the  errors  e%  are  independent,  the  joint  p.d.f.  of  the  observed  data  di  is  determined 
by  the  product  of  the  /[ej],  i.e. 

n 

f[di,  •  • .  ,dn;x,t,a\  =  JJ/[e*  =  di  -  T?:(x)  -  t\.  (2.28) 

2—1 
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Considered  as  a  function  of  the  unknown  problem  parameters,  the  data  p.d.f.  serves  as  a  like¬ 
lihood  function.  Equation  (2.25)  implies  that  the  negative  logarithm  of  the  likelihood  function 
is 

n  1 

A(x,  t,  a;  d)  =  n  log  K ( p )  +  V]  log  ^  +  n  log  cr  -\ - - tf(x,  t ;  d)  (2.29) 

~  pcrp 

1=1 

where  the  data-misfit  function  now  generalizes  equation  (2.8)  as 

n 

T(x,t;d)  =  ^2  \di  ~  T*(x)  -  t\P /vPi-  (2.30) 

i=l 

We  will  often  refer  to  the  function  A  simply  as  the  penalty  function  to  avoid  the  more  clumsy 
negative  log-likelihood. 

The  generalization  of  the  least-squares  solution  to  the  location  problem  is  the  maximum- 
likelihood  solution,  which  minimizes  the  penalty  function: 

A(x,t,<r;d)  =  min  A(x,  t,  cr;  d).  (2-31) 

x,t,cr 

We  have  included  a,  the  maximum-likelihood  estimate  of  cr,  as  part  of  the  solution,  replacing 
the  unbiased  estimate  considered  by  Flinn  (1965). 

In  practice,  the  minimization  in  equation  (2.31)  is  performed  over  a  restricted  range  of 
parameter  values.  Doing  so  effects  hard  parameter  constraints  which  can  represent  important 
prior  information,  e.g.  preventing  event  locations  above  the  Earth’s  surface  or  below  a  maximum 
credible  depth.  Additionally,  it  is  usually  credible  to  constrain  a  based  on  knowledge  of  the  data 
processing  procedures  for  picking  arrival  times.  In  the  rest  of  this  report,  it  will  be  understood 
that  the  minimization  of  a  penalty  function  is  subject  to  appropriate  a  priori  constraints  on 
the  free  parameters. 

2.3.2  Ney man- Pearson  confidence  regions 

Our  more  general  formulation  follows  the  Gaussian/linear  uncertainty  analysis  described  above 
and  defines  confidence  regions  in  terms  of  hypothesis  testing,  as  per  equations  (2.16)  and  (2.17). 
We  now  define  the  test  statistic  more  generally  as  the  logarithm  of  a  likelihood  ratio.  For  a 
confidence  region  on  x  we  have 

r(x,  d)=  min  A(x,  t,  <x;  d)  —  A(x,  f,  cr;  d),  (2.32) 

t,a 

where  A  is  the  negative  log-likelihood  function  of  equation  (2.29).  We  can  also  write  this  as 

r(x,  d)  =  A(x;  d)  —  A(x;  d),  (2.33) 

where  A  is  a  “reduced”  penalty  function  that  is,  for  any  fixed  x,  minimized  with  respect  to  t 
and  cr. 

A(x;  d)  =  minA(x,  t,  <r;d).  (2.34) 

t,cr 

We  note  that  this  statistic,  in  the  Gaussian  case  (p  =  2),  does  not  default  to  the  statistic  used 
by  Flinn  (1965)  (equation  (2.14)).  However,  the  statistics  are  equivalent  in  this  situation  and 
lead  to  the  same  confidence  region  for  any  f3. 


When  the  pick  errors  are  not  Gaussian  ( p  ^  2),  and  accepting  the  travel-time  functions 
Tj  as  nonlinear,  we  cannot  dismiss  the  possibility  that  the  distribution  of  r(x,  d)  depends  on 
some  or  all  of  the  problem  unknowns.  Therefore,  we  now  allow  rp  to  depend  on  the  unknown 
parameters  and  replace  equation  (2.16)  with 

X0(d)  =  {x:  r(x,d)  <  r^x)}  ,  (2.35) 

where 


Tfl(x)  =  maxTfl(x,t,ff).  (2.36) 

t,<7 

This  differs  from  equation  (2.16)  in  two  regards.  First,  the  dependence  of  rp  on  x  lends  a  more 
complex  nature  to  the  inequality  that  delimits  a  confidence  region.  Second,  the  dependence  on  t 
and  a  implies  that  the  null  hypothesis  being  tested  (that  x  is  the  true  hypocenter)  depends  on  t 
and  a  as  nuisance  parameters.  Equations  (2.35)  and  (2.36)  define  a  worst  case  confidence  region 
for  this  situation,  i.e.,  the  region  is  the  largest  that  any  values  for  the  nuisance  parameters  can 
yield.  A  consequence  of  this  definition  is  that  the  confidence  property  of  equation  (2.4)  is  not 
guaranteed  to  hold.  However,  a  weaker  version  of  the  property  does  hold: 


Prob 


x  €  Xp(d) 


>P- 


(2.37) 


In  this  case  one  can  state  that  the  confidence  region  Xg(d)  contains  the  true  hypocenter  with 
a  probability  of  at  least  (3. 


2.3.3  Confidence  regions  on  epicenter  and  focal  depth 

While  we  have  thus  far  presented  our  methodology  with  respect  to  hypocentral  confidence  re¬ 
gions,  and  will  continue  to  do  so,  we  pause  to  show  how  the  methodology  applies  to  confidence 
regions  on  other  event  parameters.  To  do  this,  we  separate  the  problem  parameters  into  two 
types.  The  first  are  target  parameters  whose  uncertainty  we  wish  to  characterize.  The  re¬ 
maining  parameters  are  nuisance  parameters,  which  are  not  of  interest  but  whose  effect  on  the 
data  cannot  be  ignored.  Denoting  the  target  parameters  as  components  of  a  vector  p,  and  the 
nuisance  parameters  as  components  of  q,  our  discussion  so  far  has  focused  on  the  situation 

p  =  x  (2.38) 

q  =(*,*)•  (2.39) 

Denoting  x  =  ( 9 ,  (j>,  z ),  where  9,  <j>  and  z  are  latitude,  longitude  and  depth,  respectively,  we  can 
consider  confidence  regions  on  the  event  epicenter  (9,cf>)  by  letting 

P  =  (0,<t>) 
q  =  (z,t,a), 

and  on  the  event  focal  depth  by  letting 

p  =  z  (2.42) 

q  =  (0,<t>,t,c).  (2.43) 

Our  uncertainty  formulation  transforms  to  p  and  q  as  follows.  First,  we  re-denote  the 
penalty  function  as  A(p,q;d)  and  the  test  statistic  for  hypothesis  testing  as  r(p,d).  We  have 

r(p,d)  =  A(p;d)  -  A(p;d),  (2.44) 


(2.40) 

(2.41) 
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where  (p,  q)  is  the  maximum-likelihood  solution  and 

A(p;  d)  =  min  A(p,  q;  d).  (2-45) 

p 

The  confidence  region  on  p  at  confidence  level  (3  is  then  given  by 

Vpi d)  =  {P:  T(p,d)  <  r^p)}  (2.46) 

with 

ra(p)  =  maxr/3(p,q).  (2.47) 

q 

We  will  resume  the  description  of  the  methodology  for  hypocentral  confidence  regions,  cor¬ 
responding  to  the  case  of  equations  (2.38)  and  (2.39). 

2.4  Numerical  Algorithm  for  Confidence  Regions 

Our  general  assumptions  do  not  permit  an  analytic  solution  for  tq  in  equation  (2.36)  or  a  simple 
geometric  description  of  the  confidence  region  given  by  equation  (2.35).  To  work  around  this, 
Rodi  and  Toksoz  (2000)  devised  a  numerical  approach  to  computing  and  representing  event 
location  confidence  regions,  which  we  now  describe. 

The  algorithm  involves  two  steps:  (1)  mapping  r(x,  d)  on  a  grid  in  x-space,  and  (2)  for 
given  confidence  level  (3,  estimating  7)3 (x)  for  each  point  on  the  grid.  After  performing  these 
two  steps,  the  confidence  region  at  the  (3  level  is  represented  by  flagging  the  grid  points  for 
which  t(x,  d)  <  rg(x). 

The  first  step,  which  we  will  call  likelihood  mapping,  entails  minimizing  the  penalty  function 
with  respect  to  t  and  a  with  x  fixed,  repeating  for  each  x  on  the  hypocenter  grid. 

In  this  project  we  considered  two  ways  of  performing  the  second  step,  which  we  now  discuss. 

2.4.1  Monte  Carlo  simulation 

For  any  assumed  value  of  a  (within  prescribed  bounds,  if  any),  one  can  generate  pseudo-random 
realizations  of  the  error  vector  e  based  on  the  probability  distribution  for  errors  prescribed  by 
equations  (2.25)  and  (2.27).  We  will  denote  such  a  realization  emc(<r).  Then,  for  given  x  and  t, 
a  realization  of  data  is  given  by 

dmc  =  T(x)  +  lt  +  emc(u).  (2.48) 

A  realization  of  r  is  obtained  from  this  data  as 

rmc(x,  dmc)  =  A(x;  dmc)  —  min  A(x;  dmc).  (2.49) 

X 

We  see  that  the  computation  of  rmc  entails  performing  event  location  (minimization  of  the 
penalty  function)  using  dmc  as  synthetic  data.  Generating  many  realizations  of  r  in  this  way, 
one  can  use  their  histogram  to  estimate  the  critical  statistic  value,  rg(x,  i,<r),  for  particular 
values  of  (3.  We  have  found  that  on  the  order  of  300  realizations  of  rmc  provide  adequate 
accuracy  for  rg. 

Strict  adherence  to  equations  (2.35)  and  (2.36),  however,  requires  performing  a  Monte  Carlo 
simulation,  not  just  for  each  grid  point  x,  but  for  a  sufficient  number  of  t  and  a  values  to  be  able 
to  find  a  maximum  of  rg  over  these  variables.  The  Monte  Carlo  approach  thus  becomes  very 
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computationally  intensive.  We  have  therefore  implemented  this  approach  with  the  major  sim¬ 
plification  that  the  simulation  be  performed  with  only  one  combination  of  the  parameters.  The 
parameter  values  used  are  the  maximum-likelihood  estimates.  Thus,  synthetic  data  realizations 
are  generated  as 

dmc  =  T(X)  +  l?+emc(5),  (2.50) 

and  the  formula  for  a  confidence  region  becomes 

Xp( d)  =  {x:  r(x,d)  <  T^(x,t)<r)}  .  (2.51) 

We  point  out  that  other  investigators,  notably  Wilcock  and  Toorney  (1991),  have  used  the  same 
simplification  in  similar  applications. 

Ignoring  the  dependence  of  rp  on  origin  time  is  quite  justified  since  this  dependence  arises 
only  if  bounds  on  t  are  invoked,  which  is  not  normally  done.  The  dependence  on  x  and  a  is  not 
so  easily  dismissed.  We  can  infer  from  the  Gaussian/linear  analysis  that  the  dependence  on  x 
and  a  will  be  significant  only  if  the  nonlinearity  of  the  travel-time  functions,  T),  is  significant 
or  when  the  generalize  Gaussian  order  p  differs  significantly  from  2. 

2.4.2  Likelihood  integration 

The  Bayesian  approach  to  uncertainty  suggests  an  alternative  approximation  to  773  (x).  Let 
us  assume  a  vacuous  prior  distribution  on  the  unknown  parameters,  or  /[x,  t.  a\  oc  1.  The 
likelihood  function  then  becomes  an  unnormalized  posterior  p.d.f.  on  the  parameters,  or 

f  [x,  f ,  o'  |  d]  =  JL(d)  e— A(x,t,cr’d) ,  (2.52) 

where  iv(d)  is  a  normalizing  factor.  A  Bayesian  confidence  region  on  x  (or  other  target  pa¬ 
rameter)  is  determined  by  the  marginal  posterior  p.d.f.  of  x,  as  we  indicated  earlier  in  equation 
(2.7).  Let  us  instead  consider  a  gwasi-Bayesian  approach  in  which  the  marginal  posterior  of  x 
is  replaced  with  a  distribution  that  is  maximized  with  respect  to  nuisance  parameters.  Define 

/[x  |  d]  =  max/[x,  t,  a  |  d],  (2.53) 

t,cr 

A  quasi-Bayesian  confidence  region  on  x  satisfies,  in  analogy  with  (2.7), 

/  dx  /[x  |  d]  =  fj.  (2.54) 

J  xG.Vgjd) 

Given  (2.33)  and  (2.34),  it  is  not  hard  to  show  that 

7[x  |  d]  =  K (d)  e“X(x;d)  =  K'(&)  e~T^d) ,  (2.55) 

where  K'{ d)  is  a  different  normalizing  factor.  If  we  require  that  A’g(x)  discriminate  on  the 
p.d.f.  value — i.e. 

d7(d)  =  {x:  7[x  |  d]  >  fp |  ,  (2.56) 

for  some  threshold  fp — equation  (2.55)  implies 

47(d)  =  {x :  t(x,  d)  <  t,3}  ,  (2.57) 
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True  Depth  (km)  True  Depth  (km) 


Figure  2.1:  Confidence  intervals  on  event  depth  in  a  fictitious  problem  involving  a  direct  observation  of 
depth,  which  is  either  Gaussian  (left  panel)  or  Laplace  (right)  distributed,  with  a  variance  of  (10  km)2 
in  each  case.  Three  types  of  confidence  limits  (at  (3  =  90%)  are  shown,  each  as  a  function  of  the  depth 
observation  itself  (vertical  axis):  exact  Neyman-Pearson  (green  line),  approximate  Neyman-Pearson  (red 
dots),  and  Bayesian  (blue  dots).  For  a  given  value  of  observed  depth,  the  confidence  interval  of  each 
type  consists  of  the  range  of  true  depth  between  the  confidence  limits. 


where  the  cut-off  Tg  satisfies 

[  dx  e-r(x'd)  =  p  [  dxe-r(x'd).  (2.58) 

v/r(x,d)<r(g  */t(x,  d)<oo 

This  definition  of  Tg  offers  an  alternative  approximation  to  rg  (x)  of  the  exact  Neyman-Pearson 
approach. 

A  numerical  algorithm  for  computing  quasi-Bayesian  confidence  regions  shares  the  likelihood 
mapping  step  with  the  approximate  Neyman-Pearson  method,  but  replaces  the  Monte-Carlo 
simulation  step  with  a  much  simpler  and  more  efficient  process  we  call  likelihood  integration. 
This  process  creates  a  histogram  of  the  likelihood  function  values  found  in  the  mapping  step. 
The  histogram  allows  one  to  find  approximate  solutions  of  equation  (2.58)  for  Tg  as  a  function 
of  f3.  While  the  likelihood  integration  process  is  much  faster  computationally  than  Monte  Carlo 
simulation,  it  does  place  some  additional  burden  on  the  mapping  step.  The  mapping  must  be 
sufficiently  thorough  to  yield  accurate  approximations  to  the  integrals  in  equation  (2.58). 

We  make  two  more  observations  about  approximate  Neyman-Pearson  and  quasi-Bayesian 
confidence  regions.  First,  neither  method  guarantees  the  confidence  property,  or  even  its  weaker 
version  in  equation  (2.37).  Second,  in  the  Gaussian/linear  problem,  the  methods  coincide  with 
each  other  and  with  the  exact  Neyman-Pearson  result. 

2.4.3  Comparison  in  a  simple  situation 

The  difference  between  exact  Neyman-Pearson,  approximate  Neyman-Pearson  and  quasi-Bayesian 
confidence  regions  in  the  presence  of  nonlinearity  is  illustrated  in  Figure  2.1.  This  is  a  “toy” 
problem  involving  one  unknown  parameter — event  depth — and  a  single  datum  that  is  taken  to 
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be  a  direct,  unbiased  estimate  of  the  depth.  We  can  state  this  inverse  problem  as 

d  =  z  +  e,  (2.59) 

where  d  represents  the  observed  depth.  We  introduce  nonlinearity  into  the  problem  in  the 
form  of  the  hard  constraint  that  the  true  event  depth  (horizontal  axes  in  the  figure)  cannot  be 
negative,  even  though  the  observed  depth  (vertical  axes)  might  be,  such  as  when  arrival  times 
yield  an  “airquake”  as  the  best-fitting  solution  for  an  event  hypocenter.  For  any  probability 
distribution  of  e  that  is  symmetric  about  zero  (/[e]  =  /[— e] ) ,  the  maximum-likelihood  estimate 
of  z,  subject  to  the  positivity  constraint,  is 

£  =  m  ax  (ri,  0).  (2.60) 

Intuitively,  the  positivity  constraint  should  affect  a  confidence  interval  on  depth  when  the 
observed  depth  is  negative  or  near  z  =  0,  and  we  can  see  from  the  figure  that  this  is  the  case. 
The  results  in  the  left  panel  assumes  that  the  error  in  the  observed  depth  (e)  has  a  Gaussian 
distribution  (p  =  2),  while  the  right  panel  assumes  that  e  has  a  Laplace  distribution  (p  =  1). 

The  green  lines  in  each  panel  of  Figure  2.1  trace  the  confidence  limits  for  exact  Neyman- 
Pearson  (N-P)  confidence  intervals,  which  possess  the  confidence  property.  The  red  dots  are 
the  approximate  Neyman-Pearson  confidence  limits  that  result  when  the  critical  value  Tg  is 
considered  only  for  the  maximum-likelihood  estimate  of  depth  (as  per  equation  (2.51)).  The 
blue  dots  mark  the  quasi-Bayesian  confidence  limits,  which  are  actually  proper  Bayesian  limits 
in  this  toy  problem  since  no  nuisance  parameters  are  involved.  Generally,  the  approximate 
Neyman-Pearson  intervals  are  smaller  than  the  Bayesian  intervals  for  shallow  and  negative 
depth  observations.  Which  one  better  represents  the  exact  N-P  intervals  (green  lines)  depends 
on  the  type  of  error  distribution  used. 

2.5  Examples  for  a  Nevada  Test  Site  Event 

Figure  2.2  shows  an  example  of  our  numerical  confidence  region  algorithm  applied  to  a  Nevada 
Test  Site  (NTS)  explosion  at  the  Pahute  Mesa  testing  area.  The  calculations  were  done  using 
six  Pn  arrivals  for  the  event,  assuming  Gaussian  pick  errors  with  know  variance.  The  left  panel 
shows  the  event-station  geometry.  The  center  panel  shows  the  log-likelihood  function  mapped 
on  an  epicenter  grid.  The  right  panel  shows  the  confidence  regions,  at  three  confidence  levels, 
that  result  after  performing  the  likelihood  integration  step  to  find  Tg.  The  confidence  regions 
were  computed  with  the  assumption  that  the  event  focal  depth  is  known  (fixed).  We  see  in 
the  right  panel  that  the  epicentral  confidence  region  for  each  confidence  level  is  very  close  to 
elliptical  in  shape,  indicating  that  nonlinearity  of  the  travel-time  forward  model  does  not  have 
a  significant  effect. 

Figure  2.3  shows  numerical  confidence  regions  for  the  same  Pahute  Mesa  event  computed 
with  various  non-Gaussian  pick  error  distributions.  The  confidence  regions  become  less  elliptical 
as  the  error  distribution  becomes  increasingly  less  Gaussian  (left  to  right). 
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Longitude  East  (km)  East  (km) 


Figure  2.2:  Log-likelihood  function  and  confidence  regions  for  a  Pahute  Mesa  explosion,  derived  from 
6  Pn  arrivals  with  the  pick  error  distribution  assumed  to  be  Gaussian  (p  =  2).  Left:  Schematic  map 
showing  the  event  epicenter  (red  dot)  and  station  locations  (blue  dots).  Center:  Log-likelihood  function 
mapped  on  an  epicenter  grid.  Right:  Epicenter  confidence  regions  (determined  from  the  likelihood 
function)  for  90,  95  and  98%  confidence  (blue,  green  and  red,  respectively).  In  the  right  two  panels,  the 
black  circle  marks  the  maximum-likelihood  estimate  for  the  event  epicenter,  and  the  white  circle  is  its 
GTO  epicenter  (from  Walter  et  ai,  2003). 


East  (km)  East  (km)  East  (km) 


Figure  2.3:  Numerical  confidence  regions  for  the  same  Pahute  Mesa  explosion  as  in  Figure  2.2,  computed 
with  non-Gaussian  error  distributions:  p  =  1.5  (left),  p  =  1.25  (center)  and  p  =  1  (right). 
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Chapter  3 


Model  Errors  and  the  Joint 
Location/Calibration  Inverse 
Problem 

3.1  Model  Errors 

We  extend  equation  (2.1)  to  include  an  additional  term,  ct.  which  represents  a  correction  to  the 
travel-time  predicted  by  the  function  Tt: 

di  =  Tj(x)  +t  +  Ci  +  ei,  i  =  1, . . .  ,n.  (3.1) 

Assuming  they  are  not  known,  the  travel-time  corrections  c*  play  the  role  of  prediction  errors, 
or  model  errors,  in  the  event  location  problem.  Defining  a  correction  vector, 

c  =  (ci  c2  • •  •  Cn)  T  ,  (3.2) 

we  can  write  (3.1)  in  vector  notation  as 

d  =  T(x)  +  It  +  c  +  e.  (3.3) 

Equation  (3.1)  is  a  hopelessly  ill-posed  inverse  problem  unless  we  impose  some  sort  of  prior 
constraints  on  the  Cj.  A  common  way  of  doing  this  is  to  assume  the  model  and  pick  errors 
are  both  Gaussian  random  variables  and  combine  them  as  a  single  error.  Each  data  variance, 
of,  is  then  reinterpreted  as  the  variance  of  the  error  sum,  Cj  +  e*.  This  approach  ignores  the 
different  natures  of  pick  and  model  errors  and  thus  must  be  considered  ad  hoc.  For  example, 
there  is  no  physical  basis  for  assuming  that  the  model  errors  for  different  station/phase  pairs 
are  necessarily  statistically  independent,  as  is  usually  assumed  for  pick  errors.  If  we  consider 
the  model  errors  as  distinct  from  the  pick  errors,  we  have  the  opportunity  to  constrain  them 
with  a  less  restrictive  probability  model. 

A  more  general  approach  is  to  retain  the  Cj  as  explicit  unknowns  in  the  inverse  problem 
and  constrain  them  with  prior  information  in  the  form  of  a  prior  probability  density  function. 
We  allow  this  p.d.f.  to  depend  on  the  target-event  hypocenter  and  denote  it  as  /[ c;x].  In  the 
maximum-likelihood  formulation,  the  prior  p.d.f.  on  c  is  included  by  augmenting  the  penalty 
function  for  the  arrival  data  with  a  prior  penalty  function  for  c,  given  by 

$c(c;x)  =  —  log/[c;x].  (3.4) 
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Thus,  the  likelihood  function  for  the  single-event  location  problem,  including  model  errors, 
becomes 

n  1 

A(x,  t,  c,  a;  d)  =  n\ogI\(p)  +  V'logr'i  +  nlog<r-| - -T(x,t,  c;  d)  +  <hc(c;x),  (3.5) 

P&p 

where  the  data-misfit  function  is  now 

n 

T(x,  f,  c;d)=  | dj  -  T^(x)  -  t  -  Cj\p  /uf.  (3.6) 

i= 1 

Given  this,  the  theory  of  Section  2.3  applies  to  the  new  problem  almost  as  written,  but  with 
c  adding  to  the  list  of  nuisance  parameters.  The  maximum-likelihood  solution  now  satisfies 

A(x,  t,  c,  (t;  d)  =  min  A(x,  t,  c,  a;  d).  (3.7) 

X,t,  C,CT 

The  reduced  penalty  function  is  now  given  by 

A(x;d)=  min  A(x,  t,  c,  a;  d).  (3.8) 

t,c,cr 

The  test  statistic  for  hypocenter  confidence  regions  is  still  expressed  by  equation  (2.33),  and  the 
numerical  methods  described  in  Section  2.4  still  apply.  In  the  case  of  Monte  Carlo  simulation, 
synthetic  data  are  computed  as 

dmc  =  T(X)  +  l?+c  +  emc(a).  (3.9) 

The  difficult  issue,  which  was  the  focus  of  this  project,  is  how  to  choose  the  prior  penalty 
function  <f>c(c;x).  The  approach  we  took  was  to  consider  <f>c  to  be  the  posterior  negative  log- 
likelihood  that  derives  from  a  calibration  analysis. 

3.2  Calibration  and  the  Prior  Likelihood  on  Corrections 

We  assume  that  calibration  is  performed  with  seismic  arrival-time  data  observed  from  a  set 
of  m  calibration  events.  In  analogy  with  equation  (3.1),  then,  the  calibration  problem  can  be 
expressed  as 

dij  =  Tij{xj)  +tj  +  Cij  +  eij,  i  =  1, . . .  ,rij,  j  =  1, . . .  ,m.  (3.10) 

We  do  not  assume  that  all  observable  station/phase  combinations  have  actually  been  observed 
for  every  calibration  event,  or  for  the  target  event.  Therefore,  the  particular  station/phase 
indexed  by  a  given  value  of  i  may  differ  between  events.  In  vector  notation,  equation  (3.10) 
becomes 

dj  =  T j(xj)  +  1  tj  +  cj  +  ej,  j  =  1, . . . ,  m,  (3.11) 

where,  for  a  given  j,  the  vectors  have  rij  components. 

We  also  do  not  assume  that  the  calibration  events  are  necessarily  GTO  events  having  per¬ 
fectly  known  hypocenters  and  origin  times.  In  general,  therefore,  the  unknowns  in  the  calibra¬ 
tion  inverse  problem  are  the  m  event  locations,  (xj ,  tj),  and  Y)JLi  nj  path  travel-tine  corrections, 
Cij.  With  both  correction  and  location  parameters  both  unknown,  the  calibration  problem  may 
also  be  referred  to  as  the  joint  location/calibration  problem,  although  we  reserve  that  term  for 
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the  situation  in  which  a  target  event  location  is  involved.  Additionally,  one  may  refer  to  the 
problem  as  multiple- event  location. 

As  in  the  single-event  location  problem  for  a  target  event,  we  assume  that  the  observational 
errors  in  the  calibration  problem,  e^,  have  a  generalized  Gaussian  distribution  of  order  p.  Unlike 
the  target-event  problem,  however,  we  take  the  standard  errors,  a ij,  to  be  known  completely. 
The  data  misfit  function  for  the  jth.  calibration  events  is  thus  given  by 

Uj 

(xj ,  ij,  Cj ;  dj)  =  ^2  | dij  -  Ti:j (x)  -  tj  -  Cij  |p  / (jfj .  (3. 12) 

i= 1 

3.2.1  Parameterization  of  corrections 

To  link  the  target-event  location  problem  of  equation  (3.1)  and  calibration  problem  of  equation 
(3.10),  we  must  represent  the  path  travel-time  corrections,  cl  and  c^,  with  a  common  parame¬ 
terization.  Let  us  consider  four  such  parameterizations  in  their  order  of  increasing  complexity. 
To  facilitate  this  discussion,  we  will  refer  to  the  target  event  with  the  index  j  =  0,  allowing 
equation  (3.1)  to  be  expressed  by  equation  (3.10).  Therefore,  the  corrections  c*  are  alternatively 
denoted  as  c*o,  and  the  target-event  location  parameters  (x,  t)  as  (xo,to)- 

First,  in  the  “basic”  multiple-event  location  problem  (e.g.  Jordan  and  Sverdrup,  1981;  Pavlis 
and  Booker,  1983),  travel-time  corrections  are  assumed  to  be  common  to  the  events  but  different 
between  station/phase  combinations.  If  there  are  i  unique  combinations,  the  relevant  correction 
parameters  are  time  terms,  a&,  k  =  1,  . . . ,  i,  such  that 

Cij  —  Ofcij  >  (3.13) 

where  kij  denotes  the  station/phase  pair  associated  with  the  ith  observation  for  the  jth  event. 

A  second  possible  parameterization,  introduced  by  Schultz  et  al.  (1998),  comprises  a  cor¬ 
rection  “surface,”  or  function,  specific  to  each  station/phase.  Each  surface,  a/ c(x),  defines  a 
continuous  dependence  of  the  travel-time  correction  for  an  event  on  the  event  hypocenter.  In 
this  case  we  have 

Cij  =akij(xj).  (3.14) 

A  third  example  of  a  correction  parameterization  is  the  vector  “nrislocation  function,”  a(x), 
introduced  by  Rodi  et  al.  (2005)  and  Murphy  et  al.  (2005).  In  this  case, 

Cij  =  p; ji  Xj.yij)  •  a (xj)  +  c\jj i  x , .  y , )  •  a(y^),  (3.15) 

where  y,y  is  the  station  position  associated  with  the  ijth  observation,  and  where  p,y  and  q,y  are 
slowness  vectors  (travel-time  gradients)  at  the  event  and  station  locations,  respectively.  Unlike 
equations  (3.13)  and  (3.14),  this  parameterization  ensures  source-receiver  reciprocity  (cjj  is  the 
same  if  x^  and  y,y  are  swapped). 

Finally,  we  can  relate  travel-time  corrections  to  3-D  velocity  anomalies  in  the  Earth.  Under 
the  linear  approximation  to  the  velocity  dependence  of  travel  times,  the  tomographic  parame¬ 
terization  of  travel-time  corrections  is  expressed  as 

=  J  dx  bij  (x;  Xj)  Su(x) ,  (3.16) 

where  bt]  (x;  x^)  is  a  travel-time  sensitivity  kernel  and  5u(x )  is  the  slowness  difference  between 
the  real  Earth  and  the  reference  Earth  model  used  for  evaluating  the  travel-time  function  T,j. 
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In  this  case,  the  slowness  perturbation  Su(x)  serves  as  the  common  parameterization  for  travel¬ 
time  corrections. 

Presuming  that  continuous  parameter  functions — such  as  a/ c(x),  a(x)  and  5u(x)  above — are 
discretized  on  a  spatial  grid,  we  can  generalize  these  specific  parameterizations  and  others  as 

Cij  =  bij(xj)Tu,  i  =  l,... ,  rij,  j  =  0,1, ,  m,  (3.17) 

where  the  vector  u  contains  a  discrete  set  of  parameters  for  generating  travel-time  corrections. 
It  will  be  convenient  to  write  equation  (3.17)  alternatively  as 


cj  =  Bj(xj)u,  j  =  0, 1,...  ,m, 


where  each  sensitivity  matrix  B?  (x)  is  given  by 


Bi(Xj) 


/  bij(xi)T  \ 

Mxi)T 


For  j  =  0  (target  event)  we  also  write  (3.18)  as 


(3.18) 


(3.19) 


c  =  B(x)u. 


(3.20) 


The  unknowns  in  the  calibration  inverse  problem  can  now  be  identified  as  the  parameter 
vector  u  and  calibration  event  locations  (xj,  tj),  j  =  1,  . . . ,  rn. 


3.2.2  Penalty  function  for  the  calibration  problem 

A  negative  log-likelihood  function  for  the  calibration  problem  is  given  by 


Acai(u,xi,fi, . . .  ,xm,tm;di, . . .  ,dm)  =  const 

m 

•  :  X!  x '  >  >  B/  (xi ) u;  di ) 

3= 1 
m 

+  +  <S>u(u).  (3.21) 

3= 1 


The  constant  term  (“const”)  is  not  elaborated  because  it  does  not  depend  on  any  of  the  problem 
unknowns.  The  first  summation  includes  the  data-misht  functions  for  the  calibration  arrival 
data.  These  are  given  in  equation  (3.12),  with  equation  (3.18)  substituting  for  c j.  Each  term 
of  the  second  summation  represents  a  prior  penalty  function  for  the  location  parameters  of  a 
calibration  event.  The  final  term  is  a  prior  penalty  function  for  the  correction  parameters.  An 
appropriate  choice  for  would  depend  on  the  physical  nature  of  the  parameters  in  u. 

For  this  project  we  assumed  that  prior  probability  distributions  on  x  and  t  are  of  the 
generalized  Gaussian  type.  Further,  the  epicentral  parameters  (latitude  9  and  longitude  </>), 
focal  depth  z  and  origin  time  were  assumed  to  be  independent.  Therefore,  the  prior  penalty 
function  for  a  calibration  event  is  given  by  (omitting  the  subscript  j) 


1 

const  -| - 

A($,  (f i;  9 prj ,  <f> pri ) 

PA  I 

H - 

Z  Zpri 

Pz  ! 

H - 

t  ^pri 

PA 

CTA 

Pz 

Oz 

Pt 

(3.22) 
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where  (Qpr\,  </>prj,  zpr\)  is  a  prior  hypocenter  and  tpr\  a  prior  origin  time  for  the  event.  The  function 
A  computes  the  spherical  distance  between  two  geographic  points.  Thus,  a  a  is  the  standard 
error  on  the  spherical  distance  of  the  event  epicenter  (9,(j>)  from  (0pr \,  </>prj),  while  az  and  at  are 
standard  errors  on  depth  and  time.  We  note  that  the  prior  penalty  function  effects  hard  bounds 
on  the  parameters  as  the  generalized  Gaussian  orders  (pa,  Pz,  Pt)  tend  toward  infinity. 

The  maximum-likelihood  solution  of  the  calibration  inverse  problem  is  defined  by  minimizing 
the  penalty  function  Aca|  jointly  with  respect  to  the  correction  parameters,  u,  and  the  location 
parameters  (xj,  tj)  of  the  (non-GTO)  calibration  events.  How  this  is  best  done  algorithmically 
depends  on  the  type  of  parameters  in  u.  In  designing  a  minimization  algorithm,  it  is  useful 
to  recognize,  by  examining  the  expression  for  Aca|  in  (3.21),  that,  for  any  fixed  value  of  u, 
minimization  with  respect  to  the  event  parameters  decouples  across  events  and  equates  to 
performing  single-event  location  independently  on  each  of  the  calibration  events. 


3.2.3  The  posterior  on  c 

A  prior  penalty  function  on  the  travel-time  corrections  c,  needed  for  locating  the  target  event,  is 
provided  by  a  posterior  penalty  function  on  c  that  results  from  solving  the  calibration  problem. 
In  the  maximum-likelihood  framework,  we  can  define  the  latter  in  a  quasi-Bayesian  sense, 
whereby  unknown  parameters  other  than  c  are  projected  from  the  problem  by  minimization. 
However,  since  c  itself  is  not  an  explicit  parameter  in  the  calibration  problem,  but  rather  is 
a  mapping  of  u,  projection  of  other  parameters  is  accomplished  by  constrained  minimization. 
Thus,  we  can  state  formally 


$c(c;x)=  min  min  Aca!(u, xi, H, . . .  . . ,  ,dm).  (3.23) 

u:B(x)u=c  xi,ri,...,xm,tm 


The  minimization  on  the  right-hand-side  yields  a  solution  of  the  calibration  problem  with  u 
constrained  as  B(x)u  =  c.  Substituting  from  (3.21),  we  also  have 


<f?c(c;  x)  =  const  +  min  l  y  min  p  1'&j(x.j,  tj,  Bj(xj)u;  dj) 
u:B(x)u=c  I  ^ L  i  i.i 

J=1 


+  $j(x.j,tj)  +$u(u)j. 


(3.24) 


One  way  to  interpret  this  equation  is  that  the  minimization  over  the  calibration-event  location 
parameters  transforms  the  prior  log-likelihood  on  u,  — $u(u),  into  a  posterior  log-likelihood  on 
u.  The  constrained  minimization  over  u  then  transforms  the  posterior  om  u  into  a  posterior 
on  c. 

Computing  the  function  4>c  for  use  in  locating  target  events  is  problematical.  In  general, 
<hc  does  not  admit  a  closed  form  expression.  The  alternative  of  representing  4>c  in  tabular  form 
is  not  appealing  since  a  grid  of  values  in  n-dimensional  c-space  would  be  unmanageably  large. 
Additional  difficulty  arises  from  the  dependence  of  <l?c  on  x.  The  x-dependence  also  defies  a 
closed  form  representation  and  would  add  to  the  dimensionality  of  a  tabular  representation. 


3.3  Joint  Location/Calibration  Uncertainty  Analysis 

A  confidence  region  on  the  target-event  hypocenter  x  is  obtained  by  mapping  the  function  A(x) 
in  hypocenter  space,  which  in  turn  requires  the  evaluation  of  <hc(c;  x)  for  each  grid  point  x  (see 
equations  (3.5)  and  (3.8)).  If  one  could  evaluate  equation  (3.23)  in  closed  form,  it  would  be 
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possible  to  separate  the  computations  involved  in  calibration  and  target  event  location.  In 
general,  however,  this  is  not  possible  without  invoking  some  approximation.  In  an  attempt  to 
avoid  approximations,  this  project  has  investigated  an  approach  that  bypasses  the  computation 
of  <hc  by  combining  the  target-event  location  and  calibration  problems  into  a  joint  inverse 
problem.  This  means  substituting  equation  (3.23)  into  (3.8)  and  setting  c  =  Bu,  obtaining 
(after  some  rearranging) 


A(x;  d)  =  min  <  min  A(x,  t,  B(x)u,  cr;  d) 

U  I  t,<7 


+  min  Aca|(u,xi,ti  5  •  •  •  5  -X-m  i  t"m  5  dl  5  •  •  •  ?  d.r 

XI, ti 


or 


A(x;  d)  =  const  +  min  \  min  nloga  +  p  1cr  p  min ^(x,  B(x)u;  d) 

u  L  cr  L  t 

m 

+  X!xni?  /'  '^/dx/luidp}  !  <l';ixr/;)  | 


(3.25) 


(3.26) 


The  computational  burden  resulting  from  merging  the  location  and  calibration  problems  is 
apparent.  To  evaluate  A(x;  d)  for  one  value  of  x  (i.e.  one  point  in  a  likelihood  map)  requires 
the  minimization  of  a  penalty  function  with  respect  to  all  parameters  except  x,  including  the 
calibration  parameters  (u)  and  calibration  event  locations  (xi,  t\,  ...,  xm,  tm).  That  is,  a 
large  inverse  problem  must  be  solved  repeatedly  in  order  to  calculate  a  confidence  region  on 
the  target  event  hypocenter.  The  practicality  of  the  joint  location/calibration  approach  will 
depend  on  the  nature  of  the  parameterization  of  travel-time  corrections. 


3.4  Examples  for  Basic  Multiple-Event  Location 

Computing  confidence  regions  on  a  target  event  as  part  of  a  joint  location/calibration  analysis 
is  a  feasible  task  for  the  basic  multiple-event  location  problem.  In  this  problem,  the  travel-time 
corrections  comprise  a  simple  time  term  for  each  station/phase  combination,  as  discussed  in 
Section  3.2.1  (see  equation  (3.13)).  We  have  implemented  joint  inversion  uncertainty  analysis 
for  basic  multiple-event  location  as  part  of  a  location  program,  GMEL  (Grid-search  Multiple- 
Event  Location),  developed  under  this  and  previous  projects.  The  GMEL  algorithm  is  described 
by  Rodi  et  al.  (2002)  and  Rodi  (2006). 

Figure  3.1  shows  confidence  regions  for  the  same  Pahute  Mesa  event  used  earlier  (Figures 
2.2  and  2.3),  but  now  including  the  effects  of  model  errors,  which  are  equated  to  the  uncertainty 
in  the  station/phase  time  terms  estimated  from  calibration  events.  We  used  32  other  explo¬ 
sions  at  Pahute  Mesa  and  Rainier  Mesa  as  the  calibration  events.  Only  one  of  the  calibration 
events,  a  relatively  well-recorded  event  at  Rainier  Mesa  (16  Pn  arrivals),  was  assigned  a  finite 
ground-truth  level.  The  three  panels  in  the  figure  show  the  resulting  confidence  regions  under 
three  assumptions  about  the  GT  level  of  that  calibration  event:  GTO,  GT2  or  GT5  (at  90% 
confidence).  We  see  that  the  confidence  regions,  which  now  take  model  errors  into  account,  are 
larger  than  when  model  errors  are  assumed  to  be  zero  (Figure  2.2),  and  grow  as  the  uncertainty 
in  the  GT  calibration  event  location  is  increased.  Figure  3.2  repeats  the  GTO  case  with  various 
non-Gaussian  pick  error  distribution:  p  =  1.5,  1.25  and  1  (left,  center,  right,  respectively). 

The  confidence  regions  shown  in  Figures  3.1  and  3.2  took  between  5  and  30  CPU  minutes 
each  to  compute  on  a  high-end  Linux  workstation,  the  ones  for  smaller  p  taking  the  longest.  In 
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East  (km)  East  (km)  East  (km) 


Figure  3.1:  Confidence  regions  for  the  Pahute  Mesa  target  event  shown  in  Figure  2.2,  but  accounting  for 
uncertainty  in  travel-time  corrections  (model  errors).  The  corrections  were  constrained  by  32  calibration 
events  (other  NTS  explosions)  with  one  of  them  assigned  a  finite  GT  level:  GTO  (left),  GT2  (center),  or 
GT5  (right).  The  pick  error  distribution  was  assumed  to  be  Gaussian  (p  =  2). 


East  (km)  East  (km)  East  (km) 


Figure  3.2:  Confidence  regions  for  the  same  Pahute  Mesa  event,  accounting  for  model  errors  as  con¬ 
strained  by  32  calibration  events  (see  Figure  3.1).  The  GT  calibration  event  is  assumed  to  be  GTO  and 
the  pick  error  distribution  is  non-Gaussian:  p  =  1.5  (left),  p  =  1.25  (center),  and  p  =  1  (right). 


at  least  one  case  (right  panel  in  Figure  3.2)  the  likelihood  function  was  not  mapped  sufficiently 
well  to  compute  accurate  values  of  rp,  meaning  even  more  computation  was  needed.  To  apply 
the  joint  inversion  approach  with  tomographic  corrections,  therefore,  would  be  prohibitive  since 
each  point  on  a  likelihood  grid  would  require  performing  a  full  3D  tomography  in  conjunction 
with  calibration  event  relocation. 
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Chapter  4 


Two-Stage  Uncertainty  Analysis 


To  be  computationally  practical,  in  more  than  the  simplest  cases,  the  calculation  of  calibration 
uncertainty  must  be  performed  separately  from  target-event  location.  This,  in  fact,  is  the 
traditional  practice;  the  challenge  is  to  accomplish  it  rigorously  and  efficiently.  In  the  previous 
chapter  we  mentioned  that  one  way  to  accomplish  the  separation  was  by  tabulating  a  prior 
likelihood  for  the  correction  vector  c  as  the  output  of  calibration  analysis,  and  then  using  the 
table  in  the  subsequent  location  of  a  target  event.  The  difficulty  here  is  that  the  table  would 
be  too  large;  e.g.  1014  numbers  for  20  corrections  sampled  at  five  values  each.  This  chapter 
presents  a  more  efficient  separation  scheme,  based  on  a  Gaussian  representation  of  correction 
uncertainty. 

4.1  Gaussian  (Quadratic)  Approximation 

Our  two-stage  scheme  approximates  the  prior  penalty  function  <hc(c;x)  as  quadratic  in  c.  Al¬ 
lowing  the  quadratic  function  to  be  different  for  different  x,  we  write 

<U(c;  x)  «  <hc(c*(x);  x)  +  ^(c  -  c*(x)) TVc(x)_1(c  -  c*(x)),  (4.1) 

where  Vc(x)  is  a  symmetric  matrix  and  c*(x)  is  a  stationary  point  of  <f?c: 

Vc<f>c(c  =  c*(x);  x)  =  0.  (4.2) 

When  Vc(x)  is  positive  definite — implying  4>c  is  concave  near  c  =  c*(x) — the  quadratic  ap¬ 
proximation  corresponds  to  the  assumption  that  c  has  a  Gaussian  prior  probability  distribution 
with  mean  c*(x)  and  covariance  matrix  Vc(x). 

We  now  outline  algorithms  we  developed  for  calculating  c*  and  Vc  in  a  calibration  analysis 
(Stage  1)  and  using  them  in  locating  a  target  event  (Stage  2). 

4.2  Stage  1:  Calibration 

Although  the  vector  c*  and  matrix  Vc  depend  on  x,  we  will  not  show  this  dependence  in  this 
section  for  the  sake  of  notational  clarity. 

Let  u  be  the  maximum- likelihood  solution  for  the  correction  parameter  vector  u;  i.e.,  u 
minimizes  Aca|  in  equation  (3.21)  in  conjunction  with  maximum-likelihood  estimates  of  the 
calibration  event  locations.  We  set  c*  to 

c*  =  Bu,  (4.3) 
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where  B  is  the  sensitivity  matrix  for  the  target-event  corrections.  Since  Aca|  is  a  minimum  at 
u  =  u,  equations  (3.23)  and  (4.3)  imply  that  4>c  is  a  minimum  at  c  =  c*.  In  the  absence  of 
bounds  on  c,  we  infer  that  c  =  c*  is  a  stationary  point  of  4>c. 

Our  algorithm  for  calculating  Vc  is  a  perturbation  technique.  The  algorithm  is  adapted 
from  a  method  developed  by  Rodi  and  Myers  (2007)  for  computing  travel-time  covariances 
(their  “implicit”  method);  we  refer  the  reader  to  that  paper  for  details  about  the  technique. 
Applied  here,  the  algorithm  solves  a  perturbed  calibration  problem  for  each  of  the  n  corrections, 
Cj,  associated  with  the  target  event.  Each  problem  requires  minimizing  an  augmented  penalty 
function  given  by 

A2  (u,  . . . )  =  Aca|  (u, . . . )  T  — —  (c*  +  e  n«  -  Bu) T  (c*  +  e  n«  -  Bu) .  (4.4) 

Zp 

The  vector  is  the  ith  column  of  the  identity  matrix,  containing  unity  in  its  ith  element 
and  zero  elsewhere.  The  values  of  the  scalars  e  and  p  control  the  numerical  stability  of  the 
algorithm.  Let  u®  denote  the  solution  of  the  perturbed  calibration  problem  and  define  the 
residual  vector 

rW  =  c*  _|_  enb)  _  Bu^>.  (4-5) 

Calculating  this  vector  for  each  of  the  n  perturbed  problems  we  can  construct  a  residual  matrix 
as 


R  =  (r^  •  •  •  r(n))  .  (4-6) 

Then  Vc  is  obtained  as 

Vc  =  p(eR-1-I),  (4.7) 

where  I  is  the  n  x  n  identity  matrix. 

We  note  that  the  perturbation  method  is  exact  if  Aca|,  minimized  with  respect  to  the  (x.j,tj), 
is  a  quadratic  function  of  u,  in  which  case  <hc  is  quadratic  in  c.  In  general,  however,  the  method 
can  be  viewed  as  a  technique  for  approximating  the  Hessian  matrix  of  4>c  or  fitting  a  quadratic 
function  to  <f>c. 

4.2.1  Dependence  on  x 

The  quantities  c*  and  Vc  will  depend  on  the  target  hypocenter  x  exactly  when  B,  used  in 
equations  (4.3)  and  (4.4),  depends  on  x.  The  above  algorithms  for  c*  and  Vc  are  still  valid 
when  this  dependence  occurs,  but  would  have  to  be  applied  for  each  value  of  x  of  interest. 

4.3  Stage  2:  Location  of  a  Target  Event 

Under  the  Gaussian  approximation  of  <f>c,  the  negative  log-likelihood  function  for  the  target- 
event  location  problem,  equation  (3.5),  becomes 

A(x,  t,  c,  <r;  d)  =  const  +  n  log  cr  -+ - —  \H(x,  t,  c;  d) 

+  \{c  ~  c*(x))TVc(x)-1(c  -  c*(x)),  (4.8) 
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where  T  is  given  by  equation  (3.6).  In  general,  Vc(x)  is  a  full  (non-diagonal)  matrix.  We  have 
designed  an  algorithm  to  perform  single-event  location  with  a  penalty  function  of  this  form.  To 
date,  the  algorithm  is  only  implemented  for  the  case  in  which  c*  and  Vc  do  not  depend  on  x, 
and  we  describe  this  case  first.  This  situation  occurs  with  the  time-term  parameterization  of 
travel-time  correction  discussed  in  Section  3.2.1  (equation  (3.13)),  i.e.  the  basic  multiple-event 
location  problem. 

The  algorithm  minimizes  A  in  a  hierarchical  fashion  with  the  solution  for  c  updated  in  an 
outer  loop  and  the  solution  for  the  other  parameters  found  inside  the  loop.  Let  us  define  the 
reduced  penalty  function 

A(c;d)  =  minA(x,  t,  c,  <r;d).  (4.9) 

x,£,cr 

The  outer  loop  of  the  algorithm  minimizes  A  with  a  nonlinear  conjugate  gradients  (NLCG) 
technique  (see,  for  example,  Rodi  and  Mackie,  2001).  NLCG  is  an  iterative  method  that 
produces  a  sequence  of  estimates  Co,  Ci,  . . .  that  converges  to  the  maximum- likelihood  solution, 

c. 

The  fcth  step  of  the  NLCG  iteration  requires  the  evaluation  of  A(cfc)  and  the  gradient  of  A. 
These  quantities  are  obtained  by  performing  the  minimization  in  equation  (4.9)  with  c  fixed  to 
Cfc.  That  is,  single-event  location  of  the  target  event  is  performed  with  travel-time  corrections 
fixed  to  their  current  estimate.  Our  implementation  performs  this  task  with  the  grid-search 
even  location  program  GMEL  (e.g.  Rodi,  2006).  The  result  of  the  grid-search  minimization, 
(xfc,  tk,  <Tfc),  converges  to  the  maximum-likelihood  solution  as  the  NLCG  iteration  converges. 

4.3.1  Dependence  on  x 

The  single-event  location  algorithm  just  described  could  account  for  the  dependence  of  c*  and 
Vc  on  x  in  the  grid-search  location  step,  if  these  quantities  were  computed  on  a  sufficiently 
fine  grid  in  x-space.  A  more  practical  alternative  may  be  to  embed  the  NLCG  loop  inside  a 
higher  level  loop  that  updates  the  values  of  c*  and  Vc.  The  first  NLCG  solution  would  use 
values  evaluated  at  an  initial  hypocenter  xo-  The  quantities  could  then  be  re-evaluated  at  the 
resulting  location  solution,  x,  for  use  in  a  second  NLCG  iteration.  And  so  on.  The  efficiency 
and  convergence  properties  of  such  an  approach  are  a  topic  for  future  research. 

4.4  Examples  of  Two- Stage  Approach 

We  implemented  the  two-stage  approach  for  the  problem  of  basic  multiple-event  location.  We 
show  examples  obtained  with  the  same  NTS  calibration  events,  Pahute  Mesa  target  event,  and 
Pn  arrival-time  data  used  in  previous  examples.  The  event-station  geometry  was  shown  in 
Figure  2.2. 

4.4.1  Calibration 

Using  the  perturbation  algorithm  described  in  Section  4.2,  we  computed  the  6x6  covariance 
matrix  Vc  on  the  travel-time  corrections  involved  with  the  target  event,  based  on  data  from  32 
calibration  events.  One  of  these  calibration  events  was  assigned  a  finite  GT  level  on  its  location, 
as  described  in  Section  2.5. 

Table  4.1  lists  the  correction  standard  deviations  at  two  of  the  stations  for  various  GT 
assignments.  In  addition  to  varying  the  epicenter  error  as  before  (0,  2  or  5  km),  we  also  varied 
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Table  4.1:  Correction  Standard  Deviation  (sec)  at  Stations  MNV  and  LAC 


Station  MNV 


Station  LAC 


Or.  time 


Epicenter  error  (km) 


Epicenter  error  (km) 


Correlation  Coefficient 


Figure  4.1:  The  correlation  matrix  for  travel-time  corrections  at  six  stations  recording  the  Pahute 
Mesa  target  event.  The  three  panels  correspond  to  different  GT  levels  assigned  to  the  epicenter  of  the 
ground-truth  calibration  event  at  Rainier  Mesa.  The  origin  time  error  of  the  GT  event  was  set  to  zero. 


the  origin  time  GT  level  (0,  1  or  5  sec).  The  table  shows  that  the  uncertainty  in  the  estimated 
travel-time  corrections  depends  primarily  on  the  origin-time  error  of  the  GT  event. 

Figure  4.1  displays  the  6x6  correlation  matrix  of  the  estimated  travel-time  corrections 
for  the  three  of  nine  GT  cases  in  which  the  origin  time  error  of  the  GT  event  is  zero.  Figure 

4.2  shows  the  three  cases  in  which  the  origin  time  error  is  1  second.  We  see  that  when  the 
origin-time  error  of  the  ground-truth  event  is  not  zero,  the  estimated  corrections  are  strongly 
correlated.  Not  shown  are  the  cases  where  the  GT  origin  time  error  is  5  sec;  then  the  correlations 
are  even  closer  to  1.0  than  in  Figure  4.2. 

4.4.2  Confidence  regions  on  the  target  event  epicenter 

We  generated  numerical  confidence  regions  for  the  Pahute  Mesa  target  event  using  the  travel¬ 
time  correction  means  and  covariances  obtained  from  the  calibration  stage.  The  likelihood  maps 
were  calculated  with  the  NLCG/grid-search  location  algorithm  described  in  Section  4.3.  We 
confirmed  numerically  that  confidence  regions  on  the  epicenter  do  not  depend  on  the  origin-time 
GT  level  assumed  in  calibration,  even  though  the  covariance  matrices  did. 
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ELK  LDS  NEL  LAC  MNV  BMN  ELK  LDS  NEL  LAC  MNV  BMN  ELK  LDS  NEL  LAC  MNV  BMN 


-1  -0.5  0  0.5  1 

Correlation  Coefficient 


Figure  4.2:  Same  as  Figure  4.1  except  the  origin  time  error  of  the  GT  calibration  event  is  1  second. 

Figures  4.3  and  4.4  show  confidence  regions  obtained  with  our  new  two-stage  approach  for 
the  same  cases  shown  for  the  joint  approach  in  Figures  3.1  and  3.2.  The  differences  from  the 
joint  location/calibration  results  are  small  except  for  the  GT5  case.  This  may  reflect  the  finer 
and  more  accurate  likelihood  mapping  possible  in  the  two-stage  approach  owing  to  its  greater 
efficiency.  The  likelihood  mapping  with  the  two-stage  approach  takes  only  a  few  seconds  of 
CPU  time,  rather  than  several  minutes  as  in  the  joint  approach.  Another  possible  cause  of 
the  differences  between  Figures  3.2  and  4.4  is  that  the  Gaussian  approximation  to  the  posterior 
distribution  on  travel-time  corrections  may  be  less  accurate  when  the  GT  level  of  the  calibration 
event  is  larger. 
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GTO 


GT2 


GT5 


East  (km)  East  (km)  East  (km) 


Figure  4.3:  Numerical  confidence  regions  for  the  Pahute  Mesa  target  event,  computed  with  the  two- 
stage  algorithm.  The  results  for  different  GT  levels  assigned  to  the  ground-truth  calibration  event  are 
compared  (as  in  Figure  3.1).  Gaussian  data  errors  were  assumed. 


p  =  1.5  p  =  1.25  p  =  1 


East  (km)  East  (km)  East  (km) 


Figure  4.4:  Numerical  confidence  regions  for  the  Pahute  Mesa  target  event  obtained  under  the  assump¬ 
tion  of  non-Gaussian  data  errors.  The  order  of  the  generalized  Gaussian  distribution  is  p  =  1.5  (left), 
p  =  1.25  (center)  or  p  =  1  (right).  The  epicenter  of  the  ground-truth  calibration  event  was  taken  to  be 


GTO. 
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Chapter  5 


Conclusions  and  Recommendations 


This  project  developed  a  rigorous  theoretical  framework  for  uncertainty  analysis  in  seismic 
event  location,  linking  travel-time  prediction  errors  (model  errors)  to  uncertainty  in  travel-time 
calibration.  A  complete  uncertainty  analysis  requires  considering  the  joint  location/calibration 
inverse  problem,  but  the  implementation  of  such  an  analysis,  without  simplifying  approxima¬ 
tions,  is  practical  only  for  simple  methods  of  calibration  (e.g.  the  station  time  terms  assumed 
in  basic  multiple-event  location). 

The  assumption  that  calibration  data  determine  a  Gaussian  distribution  on  travel-time 
corrections  allows  one  to  factor  the  joint  location/calibration  problem  into  two  stages.  The  cal¬ 
ibration  stage  finds  the  Gaussian  distribution,  while  the  follow-up  stage  of  target-event  location 
uses  the  distribution  as  a  prior  on  the  corrections  affecting  the  target  event.  We  have  imple¬ 
mented  and  tested  the  two-stage  approach  for  the  basic  multiple-event  location  problem  and 
confirmed  that  it  is  significantly  more  efficient  than  the  complete  joint  inversion  approach.  The 
confidence  regions  in  Figures  4.3  and  4.4  required  at  most  1  minute  of  CPU  time  for  the  calibra¬ 
tion  stage  and  less  than  5  seconds  to  compute  each  confidence  region  in  the  location  stage.  This 
indicates  that  a  rigorous  location  uncertainty  analysis  employing  complex  parameterizations  of 
travel-time  corrections  is  a  feasible  task  using  the  two-stage  approach.  Of  particular  interest 
is  the  case  of  tomographic  corrections,  parametrized  by  3-D  velocity  models.  The  recent  work 
by  Rodi  and  Myers  (2007)  provides  a  starting  point  for  implementing  the  calibration  stage  for 
this  case. 

Further  research  is  needed  to  test  the  accuracy  of  the  Gaussian  approximation  used  in  the 
two-stage  uncertainty  analysis  The  suitability  of  this  approximation  no  doubt  depends  on  the 
pick  error  distribution,  the  number  of  data  and  station  geometry,  the  GT  location  constraints 
on  calibration  events,  and  the  degree  of  nonlinearity  in  the  travel-time  forward  problem.  An 
investigation  of  this  topic  could  include  the  exploration  of  other,  more  general  approximations 
to  the  correction  posterior,  such  as  with  copulas.  Another  important  topic  for  further  study 
is  focal  depth  uncertainty.  The  formulation  and  algorithms  developed  in  this  project  apply 
to  the  computation  of  confidence  intervals  on  focal  depth  but  our  numerical  tests  considered 
only  fixed-depth  solutions  and  epicentral  confidence  regions.  Aside  from  the  importance  of 
focal  depth  in  nuclear  monitoring,  we  can  expect  focal-depth  uncertainty  to  be  significantly 
affected  by  model  errors  and  travel-time  nonlinearity,  and  therefore  take  good  advantage  of  the 
generality  of  our  approach. 
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